Simulations_network

rm(list = ls())

1 Packages

library(factoextra)
library(FactoMineR)

library(aricode)

library(kableExtra)
library(tidyverse)

library(reshape2)
library(ggplot2)
library(ggpubr)
library(gridExtra)
library(ComplexHeatmap)
library(RColorBrewer)
library(circlize)
library(dendextend)

library(qgraph)
library(igraph)
library(corrplot)

set.seed(1992)

2 Fonctions utiles

source("Functions.R")

3 Network simulation

3.1 Settings

nb_ind =  100

set.seed(2000)
classif1 = sample(1:4, size = nb_ind, replace = TRUE)
set.seed(1000)
classif2 = sample(1:3, size = nb_ind, replace = TRUE)
set.seed(10^6)
classif3 = sample(1:5, size = nb_ind, replace = TRUE)

list_clusters = list(classif1, classif2, classif3)

round(sapply(list_clusters, function(x) sapply(list_clusters, ARI, c2 = x)), 3)

##        [,1]   [,2]   [,3]
## [1,]  1.000 -0.012  0.015
## [2,] -0.012  1.000 -0.005
## [3,]  0.015 -0.005  1.000

n_net = 3 

prob_inside = 0.80 # probability of an edge between two nodes of the same cluster
prob_outside = 0.05 # probability of an edge between two nodes that are not in the same cluster

3.2 Network simulation

list_sim = unlist(lapply(list_clusters, FUN = function(groups){
  matEdgesProb = matrix(prob_outside, nrow = max(groups), ncol = max(groups))
  diag(matEdgesProb) = prob_inside
  matEdgesProb = matEdgesProb[groups, groups]
  mean_mat =  matrix(0.1, nrow = max(groups), ncol = max(groups))
  diag(mean_mat) = (1:max(groups))
  mean_mat = mean_mat[groups, groups]
  
  lapply(1:n_net, FUN = function(indice){
  
    A =  matrix(0, ncol = length(groups), nrow = length(groups))
    set.seed(indice*1000)
    A = matrix(rbinom(length(groups)^2, 1, prob = matEdgesProb), ncol = length(groups), nrow = length(groups))    
    A[lower.tri(A, diag = FALSE)] = t(A)[lower.tri(A, diag = FALSE)]
    
    return(list(A = A))
    }) 
}), recursive = FALSE)

names(list_sim) = paste0("Net_", 1:length(list_sim))

3.3 Adjancecy matrices

list_adjacency = lapply(list_sim, FUN  = function(res) res$A)

par(mfrow = c(1,3))
corrplot(list_adjacency[[1]][order(list_clusters[[1]]), order(list_clusters[[1]])], is.corr = FALSE, method = "color", diag = FALSE, bg = 'gray85')
corrplot(list_adjacency[[4]][order(list_clusters[[2]]), order(list_clusters[[2]])], is.corr = FALSE, method = "color", diag = FALSE, bg = 'gray85')
corrplot(list_adjacency[[7]][order(list_clusters[[3]]), order(list_clusters[[3]])], is.corr = FALSE, method = "color", diag = FALSE, bg = 'gray85')

4 MFA on adjacency matrices

4.1 Shortest path and MDS

dist_list = lapply(1:length(list_adjacency), FUN = function(i) distances(graph_from_adjacency_matrix(list_adjacency[[i]], mode = "undirected"), mode = "all"))

par(mfrow = c(1,3))
corrplot(dist_list[[1]][order(list_clusters[[1]]), order(list_clusters[[1]])], is.corr = FALSE, method = "color", diag = FALSE, bg = 'gray85')
corrplot(dist_list[[4]][order(list_clusters[[2]]), order(list_clusters[[2]])], is.corr = FALSE, method = "color", diag = FALSE, bg = 'gray85')
corrplot(dist_list[[7]][order(list_clusters[[3]]), order(list_clusters[[3]])], is.corr = FALSE, method = "color", diag = FALSE, bg = 'gray85')

list_mds = lapply(dist_list, cmdscale, k = nb_ind-1)

## Warning in FUN(X[[i]], ...): only 51 of the first 99 eigenvalues are > 0

## Warning in FUN(X[[i]], ...): only 54 of the first 99 eigenvalues are > 0

## Warning in FUN(X[[i]], ...): only 52 of the first 99 eigenvalues are > 0

## Warning in FUN(X[[i]], ...): only 53 of the first 99 eigenvalues are > 0

## Warning in FUN(X[[i]], ...): only 51 of the first 99 eigenvalues are > 0

## Warning in FUN(X[[i]], ...): only 53 of the first 99 eigenvalues are > 0

## Warning in FUN(X[[i]], ...): only 52 of the first 99 eigenvalues are > 0

## Warning in FUN(X[[i]], ...): only 51 of the first 99 eigenvalues are > 0

## Warning in FUN(X[[i]], ...): only 52 of the first 99 eigenvalues are > 0

4.2 MFA

res_mfa = MFA(do.call("cbind", list_mds),
              group = unlist(lapply(list_mds, ncol)),
              name.group = names(list_sim),
              type = rep("c", length(list_mds)), ncp = Inf, graph = FALSE)

fviz_eig(res_mfa, choice = "variance", addlabels = TRUE)

4.3 Heatmap Contribution

corrplot(res_mfa$group$contrib[,1:25], method = "circle", mar = c(0, 0, 1.5, 0), bg = "black", diag = TRUE, title = "Contributions, 25 first axes", is.corr = FALSE)

4.4 Factorial Maps

gg_tab = data.frame(res_mfa$group$coord, "Data" = rownames(res_mfa$group$coord), "Classification" = rep(1:3, each = n_net))
gg_tab$Classification = factor(gg_tab$Classification)

xlim_max = max(res_mfa$group$coord) + 0.05
xlim_min = -0.1
ylim_max = max(res_mfa$group$coord) + 0.05

cols = c("1" = "#E7B800", "2" = "#00AFBB", "3" = "#2E9FDF")

gg1 = ggplot(gg_tab, aes(x = Dim.1, y = Dim.2, color = Classification)) + 
  geom_point() + geom_text(aes(label=rownames(gg_tab)), size = 5,hjust=1, vjust=1) + 
  ggtitle("Axes 1/2") + theme_bw() + xlim(c(xlim_min, xlim_max)) + ylim(c(-0.05, ylim_max)) +
  xlab(paste0("Dim 1 - ", round(res_mfa$eig[1,2], 2), "%")) + ylab(paste0("Dim 2 - ", round(res_mfa$eig[2,2], 2), "%")) +
  theme(text = element_text(size = 15), axis.title = element_text(size = 15))+ 
  scale_color_manual(values = cols)
gg2 = ggplot(gg_tab, aes(x = Dim.2, y = Dim.3, color = Classification)) + 
  geom_point() + geom_text(aes(label=rownames(gg_tab)), size = 5,hjust=1, vjust=1) + 
  ggtitle("Axes 2/3") + theme_bw()  + xlim(c(xlim_min,xlim_max)) + ylim(c(-0.05,ylim_max)) +
  xlab(paste0("Dim 2 - ", round(res_mfa$eig[2,2], 2), "%")) + ylab(paste0("Dim 3 - ", round(res_mfa$eig[3,2], 2), "%")) +
  theme(text = element_text(size = 15), axis.title = element_text(size = 15))+ 
  scale_color_manual(values = cols)
gg3 = ggplot(gg_tab, aes(x = Dim.3, y = Dim.4, color = Classification)) + 
  geom_point() + geom_text(aes(label=rownames(gg_tab)), size = 5,hjust=1, vjust=1) + 
  ggtitle("Axes 3/4") + theme_bw()  + xlim(c(xlim_min,xlim_max)) + ylim(c(-0.05,ylim_max)) +
  xlab(paste0("Dim 3 - ", round(res_mfa$eig[3,2], 2), "%")) + ylab(paste0("Dim 4 - ", round(res_mfa$eig[4,2], 2), "%")) +
  theme(text = element_text(size = 15), axis.title = element_text(size = 15))+ 
  scale_color_manual(values = cols)
gg4 = ggplot(gg_tab, aes(x = Dim.4, y = Dim.5, color = Classification)) + 
  geom_point() + geom_text(aes(label=rownames(gg_tab)), size = 5,hjust=1, vjust=1) + 
  ggtitle("Axes 4/5") +  theme_bw()  + xlim(c(xlim_min,xlim_max)) + ylim(c(-0.05,ylim_max)) +
  xlab(paste0("Dim 4 - ", round(res_mfa$eig[4,2], 2), "%")) + ylab(paste0("Dim 5 - ", round(res_mfa$eig[5,2], 2), "%")) +
  theme(text = element_text(size = 15), axis.title = element_text(size = 15))+ 
  scale_color_manual(values = cols)

# postscript(file = "sim_net_group_pos.eps", width = 12, height = 3, horizontal = FALSE, onefile = FALSE, paper = "special")
# ggarrange(gg1, gg2, gg3, gg4, ncol = 2, nrow = 2, common.legend = TRUE, legend="bottom")  
# dev.off()

# png(file = "sim_net_group_pos.png", width = 1.5*19, height = 1.5*4.5, units ="cm", res = 300) 
ggall = ggarrange(gg1, gg2, gg3, gg4,  ncol = 2, nrow = 2, common.legend = TRUE, legend="bottom")
# dev.off()
ggall

hc_samples = hclust(dist(res_mfa$group$coord), method = "ward.D2")

DTC_AFM = dynamicTreeCut::cutreeDynamic(hc_samples, minClusterSize = 1, distM = as.matrix(dist(res_mfa$group$coord)))

##  ..cutHeight not given, setting it to 1.47  ===>  99% of the (truncated) height range in dendro.
##  ..done.

k_cut = max(DTC_AFM)
groups = cutree(hc_samples, k = k_cut)

# png(file = "sim_net_hc_tables.png", width = 1.5*9.5, height = 1.5*9.5, units = "cm", res = 300)
dendGraph = fviz_dend(hc_samples, k = k_cut, cex = 1.2,  k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07", "#AA4371")[1:k_cut],
          color_labels_by_k = TRUE, ggtheme = theme_bw(), horiz = TRUE, main = "Hierarchical clustering of simulated networks",
          labels_track_height = 0.3)  + theme(text = element_text(size = 14), axis.title = element_text(size = 11))
dendGraph

# dev.off()

4.5 Heatmap RV

corrplot2(res_mfa$group$RV, method = "circle", mar = c(0, 0, 0, 0), bg = "gray85", diag = FALSE, tl.cex = 1.5, tl.col = "black", 
          # title = "RV coefficients between MDS datasets",
         is.corr = FALSE, addrect = 3, order = hc_samples, addCoef.col = "antiquewhite", coef_thresh = 0.5)

corrPlot = recordPlot()

4.6 Individuals coordinates global MFA

g1_ind = fviz_mfa_ind(res_mfa, geom = "point", col.ind = factor(classif1), legend.title = "Classif 1", title = "")

## Warning in if (col.ind %in% c("cos2", "contrib", "coord")) partial = NULL: la
## condition a une longueur > 1 et seul le premier élément est utilisé

g2_ind = fviz_mfa_ind(res_mfa, geom = "point", axes = c(2,3), col.ind = factor(classif1), legend.title = "Classif 1", title = "")

## Warning in if (col.ind %in% c("cos2", "contrib", "coord")) partial = NULL: la
## condition a une longueur > 1 et seul le premier élément est utilisé

g3_ind = fviz_mfa_ind(res_mfa, geom = "point", axes = c(3,4), col.ind = factor(classif1), legend.title = "Classif 1", title = "")

## Warning in if (col.ind %in% c("cos2", "contrib", "coord")) partial = NULL: la
## condition a une longueur > 1 et seul le premier élément est utilisé

gg_classif1 = ggarrange(g1_ind, g2_ind, g3_ind, ncol = 3, common.legend = TRUE, legend = "right") 
gg_classif1 = annotate_figure(gg_classif1, top = text_grob("Individuals factorial axes - Colored according to classification 1",  vjust = 1.2, face = "bold", size = 14)) + 
  border()

g1_ind = fviz_mfa_ind(res_mfa, geom = "point", col.ind = factor(classif2), legend.title = "Classif 2", title = "")

## Warning in if (col.ind %in% c("cos2", "contrib", "coord")) partial = NULL: la
## condition a une longueur > 1 et seul le premier élément est utilisé

g2_ind = fviz_mfa_ind(res_mfa, geom = "point", axes = c(2,3), col.ind = factor(classif2), legend.title = "Classif 2", title = "")

## Warning in if (col.ind %in% c("cos2", "contrib", "coord")) partial = NULL: la
## condition a une longueur > 1 et seul le premier élément est utilisé

g3_ind = fviz_mfa_ind(res_mfa, geom = "point", axes = c(3,4), col.ind = factor(classif2), legend.title = "Classif 2", title = "")

## Warning in if (col.ind %in% c("cos2", "contrib", "coord")) partial = NULL: la
## condition a une longueur > 1 et seul le premier élément est utilisé

gg_classif2 = ggarrange(g1_ind, g2_ind, g3_ind, ncol = 3, common.legend = TRUE, legend = "right") 
gg_classif2 = annotate_figure(gg_classif2, top = text_grob("Individuals factorial axes - Colored according to classification 2",  vjust = 1.2, face = "bold", size = 14)) + 
  border()

g1_ind = fviz_mfa_ind(res_mfa, geom = "point", col.ind = factor(classif3), legend.title = "Classif 3", title = "")

## Warning in if (col.ind %in% c("cos2", "contrib", "coord")) partial = NULL: la
## condition a une longueur > 1 et seul le premier élément est utilisé

g2_ind = fviz_mfa_ind(res_mfa, geom = "point", axes = c(2,3), col.ind = factor(classif3), legend.title = "Classif 3", title = "")

## Warning in if (col.ind %in% c("cos2", "contrib", "coord")) partial = NULL: la
## condition a une longueur > 1 et seul le premier élément est utilisé

g3_ind = fviz_mfa_ind(res_mfa, geom = "point", axes = c(3,4), col.ind = factor(classif3), legend.title = "Classif 3", title = "")

## Warning in if (col.ind %in% c("cos2", "contrib", "coord")) partial = NULL: la
## condition a une longueur > 1 et seul le premier élément est utilisé

gg_classif3 = ggarrange(g1_ind, g2_ind, g3_ind, ncol = 3, common.legend = TRUE, legend = "right") 
gg_classif3 = annotate_figure(gg_classif3, top = text_grob("Individuals factorial axes - Colored according to classification 3",  vjust = 1.2, face = "bold", size = 14)) + 
  border()

4.7 Consensus network(s)

consensus_by_groups = lapply(unique(groups), FUN = function(x){
  list_adj = list_adjacency[names(list_adjacency)%in%names(groups)[groups == x]]
  Adj = Reduce("+", list_adj)
  Adj[Adj<=(0.5*length(list_adj))] = 0
  Adj[Adj>=(0.5*length(list_adj))] = 1
  Adj
})

consensus_by_groups_MFA = consensus_by_groups

round(sapply(consensus_by_groups, FUN = function(A_e){
  sapply(list_adjacency, FUN = function(A_t) sum(abs(as.matrix(A_e) - as.matrix(A_t)))/ncol(A_e)^2)
  }),2)

##       [,1] [,2] [,3]
## Net_1 0.07 0.39 0.32
## Net_2 0.07 0.39 0.32
## Net_3 0.07 0.40 0.33
## Net_4 0.38 0.08 0.36
## Net_5 0.38 0.08 0.36
## Net_6 0.39 0.08 0.36
## Net_7 0.33 0.37 0.07
## Net_8 0.33 0.38 0.07
## Net_9 0.33 0.37 0.07

res = lapply(consensus_by_groups, FUN = function(A_e){
  sapply(list_adjacency, FUN = function(A_t){
    vect = round(compareGraphs_perso(A_e, A_t), 2)
    names(vect) = c("tpr", "fpr", "tdr")
    vect})
})
names(res) = paste("Consensus ", 1:length(consensus_by_groups))

names_above = c(rep(3, 3))
names(names_above) = c(names(res))
df = do.call("rbind", res)
tab = kable(df, format = "latex") %>% kable_styling() %>% pack_rows(index = names_above)

par(mfrow = c(1,3))
qgraph(consensus_by_groups[[1]], 
       shape = "rectangle", vsize = 2, vsize2 = 2.2,
       title = paste0("Consensus on ", paste(names(groups)[groups == 1], collapse = ", ")),
       mar = c(2,2,2,2), color = c("antiquewhite", "aliceblue", "darkorchid3", "forestgreen", "blue")[classif1], diag = FALSE, 
       labels = colnames(consensus_by_groups[[1]]), title.cex = 2,
       label.cex = 1.2, negDashed = TRUE, edge.color = "black", edge.width = 0.5, layout = "spring",
       layout.par = list(repulse.rad = nrow(consensus_by_groups[[1]])^(2.5)))

qgraph(consensus_by_groups[[2]], 
       shape = "rectangle", vsize = 2, vsize2 = 2.2,
       title = paste0("Consensus on ", paste(names(groups)[groups == 2], collapse = ", ")),
       mar = c(2,2,2,2), color =  c("antiquewhite", "aliceblue", "darkorchid3", "forestgreen", "blue")[classif2], diag = FALSE, 
       labels = colnames(consensus_by_groups[[2]]), title.cex = 2,
       label.cex = 1.2, negDashed = TRUE, edge.color = "black", edge.width = 0.5, layout = "spring",
       layout.par = list(repulse.rad = nrow(consensus_by_groups[[2]])^(2.5)))

qgraph(consensus_by_groups[[3]], 
       shape = "rectangle", vsize = 2, vsize2 = 2.2,
       title = paste0("Consensus on ", paste(names(groups)[groups == 3], collapse = ", ")),
       mar = c(2,2,2,2), color =  c("antiquewhite", "aliceblue", "darkorchid3", "forestgreen", "blue")[classif3], diag = FALSE, 
       labels = colnames(consensus_by_groups[[3]]), title.cex = 2,
       label.cex = 1.2, negDashed = TRUE, edge.color = "black", edge.width = 0.5, layout = "spring",
       layout.par = list(repulse.rad = nrow(consensus_by_groups[[3]])^(2.5)))

conGraph = recordPlot()

5 Figures (MDS+MFA)

library(ggplot2)
library(cowplot)

5.1 Clustering/ Group Maps

ggsave(filename = "Figure5.pdf",  
       plot =  ggarrange(ggdraw(dendGraph, clip = "on", xlim = c(-0.05, 1)),
                         ggdraw(ggall, clip = "on"), ncol = 2,  labels = c("A)", "B)")), 
       width = 2*15, height = 2*6.5,
       units = c("cm"))

5.2 Supplementary individual coordinates

ggsave(filename = "Supplementary_unused_networks.pdf",
       plot = ggarrange(gg_classif1, gg_classif2, gg_classif3, labels = c("A)", "B)", "C)"),  ncol = 1), 
       width = 2*12, height = 2*12, units = c("cm"))

5.3 Consensus networks

png(file = "Figure6.png", width = 2*15, height = 2*5, units = "cm", res = 300)
par(mfrow = c(1,3))
qgraph(consensus_by_groups[[1]], 
       shape = "rectangle", vsize = 2, vsize2 = 2.2,
       title = paste0("Consensus on ", paste(names(groups)[groups == 1], collapse = ", ")),
       mar = c(2,2,2,2), color = c("antiquewhite", "aliceblue", "darkorchid3", "forestgreen", "blue")[classif1], diag = FALSE, 
       labels = colnames(consensus_by_groups[[1]]), title.cex = 1.5,
       label.cex = 1.2, negDashed = TRUE, edge.color = "black", edge.width = 0.9, layout = "spring",
       layout.par = list(repulse.rad = nrow(consensus_by_groups[[1]])^(2.5)))

qgraph(consensus_by_groups[[2]], 
       shape = "rectangle", vsize = 2, vsize2 = 2.2,
       title = paste0("Consensus on ", paste(names(groups)[groups == 2], collapse = ", ")),
       mar = c(2,2,2,2), color =  c("antiquewhite", "aliceblue", "darkorchid3", "forestgreen", "blue")[classif2], diag = FALSE, 
       labels = colnames(consensus_by_groups[[2]]), title.cex = 1.5,
       label.cex = 1.2, negDashed = TRUE, edge.color = "black", edge.width = 0.9, layout = "spring",
       layout.par = list(repulse.rad = nrow(consensus_by_groups[[2]])^(2.5)))

qgraph(consensus_by_groups[[3]], 
       shape = "rectangle", vsize = 2, vsize2 = 2.2,
       title = paste0("Consensus on ", paste(names(groups)[groups == 3], collapse = ", ")),
       mar = c(2,2,2,2), color =  c("antiquewhite", "aliceblue", "darkorchid3", "forestgreen", "blue")[classif3], diag = FALSE, 
       labels = colnames(consensus_by_groups[[3]]), title.cex = 1.5,
       label.cex = 1.2, negDashed = TRUE, edge.color = "black", edge.width = 0.9, layout = "spring",
       layout.par = list(repulse.rad = nrow(consensus_by_groups[[3]])^(2.5)))
dev.off()

## png 
##   2

6 Simulations with kernels

library(mixKernel)

6.1 Transformation of shortest path matrices into kernels

names(dist_list) = names(list_sim)

list_kernel_mat = lapply(dist_list, FUN = function(mat) -double_centering(mat^2)/2)

list_kernel_mat = lapply(list_kernel_mat, FUN = function(mat){
  eigens = eigen(mat, only.values = FALSE, symmetric = TRUE)
  eigens$vectors[,eigens$values>0]%*%diag(eigens$values[eigens$values>0])%*%t(eigens$vectors[,eigens$values>0])
})

6.2 Combinaison des kernels

list_kernel_mat = lapply(list_kernel_mat, FUN = function(mat){
  kern = list(kernel = mat, kernel.fun = "shortest_path")
  class(kern) = "kernel"
  kern
})

list_kernels = list("Net_1" = list_kernel_mat$Net_1,
           "Net_2" = list_kernel_mat$Net_2,
           "Net_3" = list_kernel_mat$Net_3,
           "Net_4" = list_kernel_mat$Net_4,
           "Net_5" = list_kernel_mat$Net_5,
           "Net_6" = list_kernel_mat$Net_6,
           "Net_7" = list_kernel_mat$Net_7,
           "Net_8" = list_kernel_mat$Net_8,
           "Net_9" = list_kernel_mat$Net_9)

list_args = list_kernels
list_args$method = "full-UMKL"
combKern = do.call(combine.kernels, list_args)

round(combKern$weights, 3)

## [1] 0.111 0.111 0.111 0.111 0.111 0.111 0.111 0.111 0.111

combKern$weights

## [1] 0.1111022 0.1110528 0.1112112 0.1112315 0.1112853 0.1112747 0.1109838
## [8] 0.1109184 0.1109400

mat_sim = similarities(list_kernels)
rownames(mat_sim) = colnames(mat_sim) = names(list_kernels)

mat_dist = DistFromSim(mat_sim)

hc_samples = hclust(as.dist(mat_dist), method = "complete")

corrplot2(mat_sim, method = "circle", mar = c(0, 0, 0, 0), bg = "gray85", diag = FALSE, tl.cex = 1, tl.col = "black", 
         is.corr = FALSE, addrect = 3, order = hc_samples, addCoef.col = "antiquewhite", coef_thresh = 0.5)

corrPlot = recordPlot()

hc_samples = hclust(as.dist(DistFromSim(mat_sim)), method = "complete")

DTC_Kern = dynamicTreeCut::cutreeDynamic(hc_samples, minClusterSize = 1, distM = as.matrix(as.dist(mat_dist)))

##  ..cutHeight not given, setting it to 1.26  ===>  99% of the (truncated) height range in dendro.
##  ..done.

k_cut = max(DTC_Kern)
groups = cutree(hc_samples, k = k_cut)

# png(file = "sim_hc_net_Kernel.png", width = 1.5*9.5, height = 1.5*9.5, units = "cm", res = 300)
fviz_dend(hc_samples, k = 3, cex = 1.5,  k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07", "#AA4371"),
          color_labels_by_k = TRUE, ggtheme = theme_bw(), horiz = TRUE, main = "Hierarchical clustering of simulated graphs",
          labels_track_height = 0.5) + theme(text = element_text(size = 17), axis.title = element_text(size = 15))

## Warning in get_col(col, k): Length of color vector was longer than the number of
## clusters - first k elements are used

# dev.off()

# postscript(file = "sim_hc_tables.eps", width = 8, height = 8, horizontal = FALSE, onefile = FALSE, paper = "special")
dendGraph = fviz_dend(hc_samples, k = 3, cex = 1,  k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07", "#AA4371")[1:3],
          color_labels_by_k = TRUE, ggtheme = theme_bw(), horiz = TRUE, main = "Hierarchical clustering of simulated networks",
          labels_track_height = 0.5) + theme(text = element_text(size = 12), axis.title = element_text(size = 12)) + 
  theme(plot.margin = unit(c(0,0.1,0,0.1), "cm"), aspect.ratio = 1)
dendGraph

# dev.off()

ggsave(filename = "Figure7.pdf",  
       plot =  ggarrange(ggdraw(dendGraph, clip = "on", xlim = c(-0.05, 1)),  ggdraw(corrPlot, clip = "on"), ncol = 2, labels = c("A)", "B)")),
       width = 2*15, height = 2*6.5, units = c("cm"))

6.3 Consensus network(s)

consensus_by_groups = lapply(unique(groups), FUN = function(x){
  list_adj = list_adjacency[names(list_adjacency)%in%names(groups)[groups == x]]
  Adj = Reduce("+", list_adj)
  Adj[Adj<=(0.5*length(list_adj))] = 0
  Adj[Adj>=(0.5*length(list_adj))] = 1
  Adj
})

round(sapply(consensus_by_groups, FUN = function(A_e){
  sapply(consensus_by_groups_MFA, FUN = function(A_t) sum(abs(as.matrix(A_e) - as.matrix(A_t)))/ncol(A_e)^2)
  }),2)

##      [,1] [,2] [,3]
## [1,] 0.00 0.38 0.31
## [2,] 0.38 0.00 0.36
## [3,] 0.31 0.36 0.00

par(mfrow = c(1,3))
qgraph(consensus_by_groups[[1]], 
       shape = "rectangle", vsize = 2, vsize2 = 2.2,
       title = paste0("Consensus on ", paste(names(groups)[groups == 1], collapse = ", ")),
       mar = c(2,2,2,2), color = c("antiquewhite", "aliceblue", "darkorchid3", "forestgreen", "blue")[classif1], diag = FALSE, 
       labels = colnames(consensus_by_groups[[1]]), title.cex = 2,
       label.cex = 1.2, negDashed = TRUE, edge.color = "black", edge.width = 0.5, layout = "spring",
       layout.par = list(repulse.rad = nrow(consensus_by_groups[[1]])^(2.5)))

qgraph(consensus_by_groups[[2]], 
       shape = "rectangle", vsize = 2, vsize2 = 2.2,
       title = paste0("Consensus on ", paste(names(groups)[groups == 2], collapse = ", ")),
       mar = c(2,2,2,2), color =  c("antiquewhite", "aliceblue", "darkorchid3", "forestgreen", "blue")[classif2], diag = FALSE, 
       labels = colnames(consensus_by_groups[[2]]), title.cex = 2,
       label.cex = 1.2, negDashed = TRUE, edge.color = "black", edge.width = 0.5, layout = "spring",
       layout.par = list(repulse.rad = nrow(consensus_by_groups[[2]])^(2.5)))

qgraph(consensus_by_groups[[3]], 
       shape = "rectangle", vsize = 2, vsize2 = 2.2,
       title = paste0("Consensus on ", paste(names(groups)[groups == 3], collapse = ", ")),
       mar = c(2,2,2,2), color =  c("antiquewhite", "aliceblue", "darkorchid3", "forestgreen", "blue")[classif3], diag = FALSE, 
       labels = colnames(consensus_by_groups[[3]]), title.cex = 2,
       label.cex = 1.2, negDashed = TRUE, edge.color = "black", edge.width = 0.5, layout = "spring",
       layout.par = list(repulse.rad = nrow(consensus_by_groups[[3]])^(2.5)))

conGraph = recordPlot()

7 InfoSession

sessionInfo()

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
## 
## locale:
##  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
##  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
##  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] mixKernel_0.4        reticulate_1.16      mixOmics_6.12.2     
##  [4] lattice_0.20-41      MASS_7.3-53          cowplot_1.1.0       
##  [7] corrplot_0.84        igraph_1.2.5         qgraph_1.6.5        
## [10] dendextend_1.14.0    circlize_0.4.10      RColorBrewer_1.1-2  
## [13] ComplexHeatmap_2.4.3 gridExtra_2.3        ggpubr_0.4.0        
## [16] reshape2_1.4.4       forcats_0.5.0        stringr_1.4.0       
## [19] dplyr_1.0.2          purrr_0.3.4          readr_1.3.1         
## [22] tidyr_1.1.2          tibble_3.0.3         tidyverse_1.3.0     
## [25] kableExtra_1.2.1     aricode_1.0.0        FactoMineR_2.3      
## [28] factoextra_1.0.7     ggplot2_3.3.2       
## 
## loaded via a namespace (and not attached):
##   [1] tidyselect_1.1.0      htmlwidgets_1.5.1     munsell_0.5.0        
##   [4] codetools_0.2-16      withr_2.3.0           colorspace_1.4-1     
##   [7] Biobase_2.48.0        phyloseq_1.32.0       knitr_1.30           
##  [10] rstudioapi_0.11       leaps_3.1             stats4_4.0.3         
##  [13] ggsignif_0.6.0        huge_1.3.4.1          labeling_0.3         
##  [16] mnormt_2.0.2          farver_2.0.3          rhdf5_2.32.2         
##  [19] vctrs_0.3.4           generics_0.0.2        xfun_0.17            
##  [22] R6_2.4.1              clue_0.3-57           gridGraphics_0.5-0   
##  [25] assertthat_0.2.1      scales_1.1.1          nnet_7.3-14          
##  [28] gtable_0.3.0          rlang_0.4.7           scatterplot3d_0.3-41 
##  [31] GlobalOptions_0.1.2   splines_4.0.3         rstatix_0.6.0        
##  [34] broom_0.7.0           checkmate_2.0.0       yaml_2.2.1           
##  [37] abind_1.4-5           modelr_0.1.8          d3Network_0.5.2.1    
##  [40] backports_1.1.10      Hmisc_4.4-1           tools_4.0.3          
##  [43] psych_2.0.8           lavaan_0.6-7          ellipsis_0.3.1       
##  [46] biomformat_1.16.0     BiocGenerics_0.34.0   dynamicTreeCut_1.63-1
##  [49] Rcpp_1.0.5            plyr_1.8.6            base64enc_0.1-3      
##  [52] zlibbioc_1.34.0       rpart_4.1-15          pbapply_1.4-3        
##  [55] GetoptLong_1.0.2      viridis_0.5.1         S4Vectors_0.26.1     
##  [58] haven_2.3.1           ggrepel_0.8.2         cluster_2.1.0        
##  [61] fs_1.5.0              magrittr_1.5          data.table_1.13.0    
##  [64] RSpectra_0.16-0       openxlsx_4.2.2        reprex_0.3.0         
##  [67] tmvnsim_1.0-2         whisker_0.4           matrixStats_0.57.0   
##  [70] hms_0.5.3             evaluate_0.14         rio_0.5.16           
##  [73] jpeg_0.1-8.1          readxl_1.3.1          IRanges_2.22.2       
##  [76] shape_1.4.5           compiler_4.0.3        ellipse_0.4.2        
##  [79] crayon_1.3.4          htmltools_0.5.0       mgcv_1.8-33          
##  [82] corpcor_1.6.9         Formula_1.2-3         lubridate_1.7.9      
##  [85] DBI_1.1.0             dbplyr_1.4.4          Matrix_1.2-18        
##  [88] ade4_1.7-15           car_3.0-9             permute_0.9-5        
##  [91] cli_2.0.2             quadprog_1.5-8        parallel_4.0.3       
##  [94] BDgraph_2.63          pkgconfig_2.0.3       flashClust_1.01-2    
##  [97] foreign_0.8-80        xml2_1.3.2            foreach_1.5.0        
## [100] pbivnorm_0.6.0        rARPACK_0.11-0        multtest_2.44.0      
## [103] webshot_0.5.2         XVector_0.28.0        LDRTools_0.2-1       
## [106] rvest_0.3.6           digest_0.6.25         vegan_2.5-6          
## [109] Biostrings_2.56.0     rmarkdown_2.3         cellranger_1.1.0     
## [112] htmlTable_2.1.0       curl_4.3              gtools_3.8.2         
## [115] rjson_0.2.20          lifecycle_0.2.0       nlme_3.1-149         
## [118] glasso_1.11           jsonlite_1.7.1        Rhdf5lib_1.10.1      
## [121] carData_3.0-4         viridisLite_0.3.0     fansi_0.4.1          
## [124] pillar_1.4.6          httr_1.4.2            survival_3.2-7       
## [127] glue_1.4.2            zip_2.1.1             fdrtool_1.2.15       
## [130] png_0.1-7             iterators_1.0.12      stringi_1.5.3        
## [133] blob_1.2.1            latticeExtra_0.6-29   ape_5.4-1